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ABSTRACT 

Context. Dynamos in the Sun and other bodies tend to produce magnetic fields that possess magnetic helicity of opposite sign at 
large and small scales, respectively The build-up of magnetic helicity at small scales provides an important saturation mechanism. 
Aims. In order to understand the nature of the solar dynamo we need to understand the details of the saturation mechanism in spherical 
geometry. In particular, we want to understand the eff'ects of magnetic helicity fluxes from turbulence and meridional circulation. 
Methods. We consider a model with just radial shear confined to a thin layer (tachocline) at the bottom of the convection zone. The 
kinetic a owing to helical turbulence is assumed to be localized in a region above the convection zone. The dynamical quenching 
formalism is used to describe the build-up of mean magnetic helicity in the model, which results in a magnetic a eff'ect that feeds back 
on the kinetic a eff'ect. In some cases we compare with results obtained using a simple algebraic a quenching formula. 
Results. In agreement with earlier findings, the magnetic a eff'ect in the dynamical a quenching formalism has the opposite sign 
compared with the kinetic a eff'ect and leads to a catastrophic decrease of the saturation field strength with increasing magnetic 
Reynolds numbers. However, at high latitudes this quenching eff'ect can lead to secondary dynamo waves that propagate poleward due 
to the opposite sign of a. Magnetic helicity fluxes both from turbulent mixing and from meridional circulation alleviate catastrophic 
quenching. 
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1. Introduction 



The solar dynamo models developed so far and which agree with 
solar magnetic field observations tend to solve the aQ mean 
field dynamo equations. The turbulent a-eflTect first proposed 
by Parker (1955) is believed to be generated due to helical tur- 
bulence in the convection zone of the Sun. Since a is gener- 
ated due to quadratic correlations of the small-scale turbulence 
we need a closure in order to complete the set of mean field 
equations, e.g., the first order smoothing approximation (FOSA), 
and express the mean electromotive force in terms of the mean 
magnetic fields. This turbulent a encounters a critical problem 
when the energy of the mean field becomes comparable to the 
equipartition energy of the turbulence in the convection zone 
and hence it becomes increasingly difl^icult for the helical turbu- 
lence to twist rising blobs of magnetic field. The solar dynamo 
modellers have traditionally used what is referred to as algebraic 
alpha quenching to mimic this phenomena. This involves replac- 

ing ahy ao/(l-\-B /B^^), an expression used since Jepps (1975), 
or by ao/(l + Rj^B /B^^), where a^ is the unquenched value and 

Rj^ is the magnetic Reynolds number, B is the mean magnetic 
field and B^^ is the equipartition magnetic field. The latter ex- 
pression has been discussed since the early work of Vainshtein 
& Cattaneo (1992). The Rj^ in the denominator comes from 
the fact that the small-scale fluctuating magnetic field reaches 
equipartition long before the mean magnetic field does. This has 
been supported by several numerical experiments to determine 
the saturation behaviour of a (e.g. Cattaneo & Hughes 1996, 
Ossendrijver et al. 2002). Given the large magnetic Reynolds 
numbers of Astronomical objects, such phenomena is referred 
to as catastrophic quenching. 



After the discovery of the layer of strong radial shear (called 
the tachocline by Spiegel & Zahn 1992) at the bottom of the so- 
lar convection zone, Parker (1993) proposed a new class of solar 
dynamo models called the interface dynamo. In these models 
the shear is confined to a narrow overshoot layer just beneath 
the convection zone, also the region of a eflTect. The dynamo 
wave propagates in a direction given by the Parker- Yoshimura 
rule at the interface between the two layers defined by a steep 
gradient in the turbulent diflTusivity. The toroidal field produced 
due to stretching by the shear is much stronger than the poloidal 
field and remains confined in the overshoot layer, away from 
the region where the a eflTect operates. It may be noted that the 
interface dynamo model may have serious problems when solar- 
like rotation with positive latitudinal shear is included (Markiel 
& Thomas 1999). Similarly, in the Babcock-Leighton class of 
flux transport models (Choudhuri et al. 1995; Durney 1995) the 
toroidal and the poloidal fields are produced in two diflTerent 
layers. Unlike in the interface dynamo models, the coupling be- 
tween the two layers is mediated both by diflTusion and the con- 
veyer belt mechanism of the meridional circulation. 

It has been proposed that in interface and Babcock-Leighton 
type dynamos, the a eflTect is not catastrophically quenched at 
high Rj^ because the strength of the toroidal field is very weak in 
the region of finite turbulent a (e.g. Tobias, 1996; Charbonneau, 
2005). However, according to our knowledge, not much has been 
done to study the variation of the amplitude of the saturation 
magnetic field with the magnetic Reynolds number for these 
classes of afl dynamos. Zhang et al (2006) made an attempt 
to reproduce the surface observations of current helicity in the 
Sun using a 2D mean field dynamo model in spherical coordi- 
nates coupled with the dynamical quenching equation. In a sepa- 
rate paper (Chatterjee, Brandenburg & Guerrero, 2010) we have 



p. Chatterjee et al.: Magnetic helicity fluxes in spherical shell dynamos 



demonstrated that interface dynamo models are also subject to 
catastrophic quenching. 

It has been identified a decade ago that the small-scale mag- 
netic helicity generated due to the dynamo action back reacts on 
the helical turbulence and quenches the dynamo (Blackman & 
Field, 2000; Kleeorin et al. 2000). It has now been shown that 
this mechanism reduces the saturation amplitude of the mag- 
netic field (Bsat) with increasing magnetic Reynolds number 
(Rm)- Nevertheless this constraint may be lifted if the system 
is able to get rid of small scale helicity through several ways 
like open boundaries, advective, diff'usive and shear driven fluxes 
(Shukurov et al. 2006, Zhang et al. 2006, Sur et al. 2007, Kapyla 
et al. 2008, Brandenburg et al. 2009, Guerrero et al. 2010). Even 
though the helicity constraint in direct numerical simulations 
(DNS) of dynamos with strong shear have been clearly identi- 
fied, the results can be matched with mean field models having 
a weaker algebraic quenching than a^ dynamos (Brandenburg 
et al. 2001). It is possible to include this process in mean-field 
dynamo models through an equation describing the evolution of 
the small scale current helicity. We shall refer to this equation as 
the dynamical quenching mechanism. 

In this paper we perform a series of calculations with mean 
field QfQ models in spherical geometry along with a dynamical 
equation for the evolution of a for magnetic Reynolds numbers 
in the range 1 < 7?^ < 2 x 10^ . An important feature of the calcu- 
lation is that the region of strong narrow shear is separated from 
the region of helical turbulence. This paper in addition to provid- 
ing detailed results not mentioned in Chatterjee, Brandenburg & 
Guerrero (2010), is also aimed at studying somewhat more com- 
plicated models including meridional circulation. The role of 
diff'usive helicity fluxes modelled into the dynamical quenching 
equation by using a Fickian difl'usion term is also discussed for 
various models. It may be mentioned that helicity fluxes across 
an equator can indeed be modelled by such a difl'usion term as 
shown by Mitra et al. (2010). In §2 we discuss the features of the 
QfQ model used, and the formulation of dynamical a quenching. 
The results are highlighted in §3 and conclusions are drawn in 
§4. 

2. The basic aQ Dynamo Model 

2.1. Simple two-layer dynamo 

We solve the induction equation in a spherical shell assuming 
axisymmetry. Our dynamo equations consists of the induction 
equations for the mean poloidal potential A(p(r, 6) and the mean 
toroidal field ^^(r, 6). Axisymmetry demands that for all vari- 
ables didcp = 0. Let us first do a qualitative estimate of the turbu- 
lent a and the turbulent diff'usivity %• From mixing length theory 
we have (cf. Sur et al. 2008), 



% 



'^rms 



where u^ms is the rms velocity of the turbulent eddies, kf is the 
wavenumber of the energy-carrying eddies, corresponding to the 
inverse pressure scale height near the base of the convection 
zone. Since we have made use of the error function profile ex- 
tensively, let us denote 



0-(r,rc,4) = l±erf 



r - Tg 



We have used a smoothed step profile for rj^ given by 

T](r) = rir + r]t&^(r, r^, de) (1) 



where r^ = 0.737?©, and de = 0.0257?©. In this paper we define 
the magnetic Reynolds number Rj^ = rjt/rjr. Using FOSA we 
also have ao = T6f6L)rmsWrms/3, where oj^ms is the rms vorticity 
of the turbulence and r ~ (kfU^j^s)'^ is the eddy correlation time 
scale. The prefactor 6f, usually of order 0.1 or less is used since 
(u.co\ms < Wrms^rms- The casc 6f = 1 mcans the flow is maxi- 
mally helical. These approximations give us an estimate of ao in 
terms of eddy difl'usivity T]t and forcing scale kf as. 



ao = €f- 



rkful^ 



= ^mh. 



We would consider kf rather than a^ as a free parameter in the 
model apart from z/t. Assuming equipartition between magnetic 
energy and the turbulent energy, we also calculate an equiparti- 
tion magnetic field ^eq as, 

^eq = {^npfl^u^, = (47ip)'^^3r],kf. 

For algebraic quenching we consider the following form for 
kinematic ckk given by. 



QfK(r) = 



O.56f77t^f0^(r, r^, da) cos sin^ 6 



-^^ , 



(2) 



1 + g,B IBi^ 



where ga is a non-dimensional coefficient equal to 1 or R^ de- 
pending on the assumed form of algebraic quenching in the mod- 
els and ^ = unless given. Even though the helical turbulence 
pervades almost the entire convection zone, we take r^ = 0.777?© 
and da - 0.0157?© so that we can have a large separation between 
the shear and turbulent layer. Consequently we consider a diff'er- 
ential rotation profile like that in the high latitude tachocline of 
the Sun given by. 



Q(r) = -Do0^(r,r^,^w), 



(3) 



where Qq = 14nHz, r^ = 0.687?© and d^ = 0.0157?©. The radial 
profiles of 7/t, a and dO^ldr are plotted as a function of fractional 
radius r/7?© in Fig. [T] The region of strong radial shear is sep- 
arated from the region of helical turbulence and the diff'usivity 
has a strong gradient at a radius lying between the layers of fi- 
nite strong shear and turbulent a. The reason of the same is to 
decrease the time period T^yx of the oscillatory dynamos to a 
reasonably small fraction of the diff'usion time ^diff • Otir aim is 
to solve the induction equations coupled with yet another equa- 
tion for the evolution of a-eff'ect, the formulation of which is 
described in §2.1. 

2.2. Dynamical a quenching 

It was first shown by Pouquet et al. (1976) that the turbulent a 
eff'ect is modified due to the generation of small-scale helicity in 
the way given by Eq. ^ below. The second term is sometimes 
referred to as the magnetic of-eff'ect. 



a = (2k + C^M 



a> ' u - p ^j 



T). 



(4) 



where co, u, j\ b denote the fluctuating component of the vor- 
ticity, velocity, current and magnetic field in the plasma. It is 
possible to write an equation for the evolution of the magnetic 
part of Of or au from the equation for evolution of the small-scale 
magnetic helicity density hf = a ' b using the relation. 



au 



^eq 



hf. 



(5) 



p. Chatterjee et al.: Magnetic helicity fluxes in spherical shell dynamos 



L2 

0.8 
0,6 

0,4 b 

0,2 
0,0 



dQ/dr 



I \ 

I \ 

I \ 

I \ 

I \ 

I \ 

I \ 

I \ 

I \ 

I \ 




dial and the latitudinal components of the meridional flow, Up 
here. 



0.60 0.65 0.70 0.75 0.80 0.85 
r/Rs 

Fig. 1. Profiles of radial shear dfl/dr (nHz cm"^), a (cm s"^) 
and T] (10^^ cm^ s"^) as a function of fractional solar radius. 



However the equation for a ' b will be gauge-dependent and 
it makes sense only to write an equation for the volume aver- 
aged quantity in order to avoid dependence on specific gauge 
(Blackman & Brandenburg 2002). Our dynamo equations are 
independent of any gauge since we solve for the magnetic po- 
tential component A^ with an axisymmetric constraint. It is im- 
portant for us that the equation for au is also gauge independent. 
Subramanian & Brandenburg (2006) used the Gauss linking for- 
mula for the expression for hf and wrote an equation independent 
of the gauge for the magnetic helicity density under the assump- 
tion that the correlation length for all the fluctuating variables 
remain small compared to the system size at all times. Using 
Eq. ^ we write the same equation in terms of au. 



dam 
dt 



-2r]ik] 



S • B au 



V /^„ 



(6) 



where S and B are the mean field EMF and the mean magnetic 
field. The flux Fa consists of individual components, e.g., ad- 
vection due to the mean flow, Vishniac-Cho fluxes (Vishniac & 
Cho 2001), efl'ects of mean shear, difl'usive fluxes, etc. In this 
paper we have put F^, = unless mentioned otherwise. 

The decay time in Eq. (pi) is r^, = Rj^/rjik^ = 4.55 x 
lO'^Rmtdis- It should be noted that we use ga = in Eq. ^ 
whenever we employ the dynamical quenching equation, be- 
cause dynamical quenching is usually more important. 

2.3. Flux transport Babcock-Leighton dynamo 

Axisymmetric mean field solar dynamo models including merid- 
ional circulation and Babcock-Leighton a eff'ect have been stud- 
ied extensively by several authors (Dikpati & Charbonneau 
1999; Chatterjee et al. 2004; Guerrero & Dal Pino 2008, and ref- 
erences therein). These models have now reached a stage where 
they are able to reproduce the butterfly diagram and the correct 
phase between the polar fields and the toroidal fields. In this sec- 
tion we will use a Babcock-Leighton (BL) a along with an an- 
alytical meridional circulation (MC) which is poleward at the 
surface with a maximum amplitude of uq = 20 m s"^ and the 
expression for which is given by van Ballegooijen & Choudhuri 
(1988). For completeness we provide the expressions for the ra- 



Uf — Uq 



Ue = Uq 



[^)\i^-l + ^f^/' - ^^'^'](2cos'6-sm'6l (7) 



3 2 ^ 



(^) (-l+c,if^/2-c,2f^/4)sin^cos^. 



(8) 



1/2 



where ^ = RJr - 1, r^ = O.717?o, (b = Re/rb - h c,i = % 
and Cs2 = ^^l^^^ • It should be mentioned that, unlike in flux 
transport dynamo models, the meridional circulation does not 
reverse the direction of propagation of the dynamo wave in in- 
terface dynamo models as long as the meridional circulation is 
confined within the convection zone (Petrovay & Kerekes 2004). 
We solve this model along with the equation for dynamical a 
quenching described in Sect. 2.1. The fluxes in Eq. ^ are now 
given by. 



Fa = auUp - V • (KVau), 



(9) 



where k is the difl'usion coeflftcient for au taken to be KoT](r). It 
may be remembered that the ckk is now not due to the helical 
turbulence in the bulk of the convection zone, but due to a phe- 
nomenological BL a where the poloidal field is produced from 
the toroidal field due to decay of tilted bipolar active regions. 
The analytical expression for ckk is given by 

c^K = -aBL0^(r, 0.95/?o, d) ©"(r, Rq, d) cos 6 sin^ 6 (10) 

with d = O.O157?0. The BL a is assumed to be concentrated 
only in the upper 0.05% of the convection zone. The turbulent 
diff'usivity has the same profile as in Eq. ^ but with rjt = 2x 10^^ 
cm s"^ and r^ = 0.7 Rq. The shear is still radial and given by 
Eq. ^ with Tw = 0.1 Rq. 

Our computational domain is defined to be the region con- 
fined by < ^ < TT and 0.557?© < r < Rq. Unless otherwise 
stated, the boundary conditions for A^ are given by a poten- 
tial field condition at the surface (Dikpati & Choudhuri 1994) 
and A^ = at the poles. We have also performed some cal- 
culations with the vertical field condition at the top boundary, 
which means that Eq = B^p = 0. At the bottom we use the per- 
fect conductor boundary condition of Jouve et al. (2008) with 
A(p = d(rB(p)/dr = 0. However a more realistic perfect conduc- 
tor boundary condition in our opinion would be d(rBo)/dr = 
d(rB(p)/dr = 0. Also B^p = on all other boundaries. The equa- 
tion for au is an initial value problem for F^ = 0. For finite 
fluxes we have also set au = at all boundaries. We have 
checked that the results are not very sensitive to the difl'erent 
boundary conditions given above mainly because the boundaries 
are far removed from the dynamo region. 



3. Results 

3.1. Magnetic field properties without helicity fluxes 

In order to study the Rj^ dependence of the saturation magnetic 
field in the two layered dynamo with diff'usive coupling we keep 
all the dynamo parameters the same for all the runs and change 
rjr from 2 x 10^ cm^ s"^ to 2 x 10^^ cm^ s"^ while keeping rjt 
fixed at 4 X 10^^ cm^ s"^ It may also be noted that the time 
period of the dynamo models (rcyi) is fairly independent of the 
magnetic Reynolds number. We show the magnetic energies as 
a function of time for the nonlinear system with ^o = O.O^rj^kf 
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figure. 
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Fig. 3. Critical a in terms of a fraction of rj^kf as a function of 
magnetic Reynolds number Rj^ for the interface dynamo model 
of Fig. [1] 



for a range of magnetic Reynolds numbers in Fig. [2] The strong 
Rjn dependence which is reminiscent of catastrophic quenching 
in all astrophysical dynamos can be easily discerned from Fig.[^ 
It is interesting that the saturation energy of the Rm = ^ model is 
lower than that of the Rj^ = 20 case. The dynamo model may be 
highly dissipative at very low magnetic Reynolds numbers. 

The slopes of the volume averaged energy are also very dif- 
ferent in the kinematic phase, which means that the critical dy- 
namo numbers also depend on Rj^. To be able to correctly com- 
pare the dynamo models for diff'erent Rj^, it is first important to 
calculate the critical value of ao, denoted by ac for each model. 
Such a plot is shown in Fig. |3] From this figure we can con- 
clude that this dynamo model is most efficient near Rj^ = 20. 
A similar variation of ac with the ratio ^/t/^r was obtained an- 
alytically for interface dynamos by MacGregor & Charbonneau 
(1997; see their Fig. 5A). We now set ao = 2a c, corresponding 
to the /?ni of each model, and repeat our calculations. We shall 



now use this value of a for the rest of the paper. The satura- 
tion energy decreases monotonically as a function of magnetic 
Reynolds number as shown in Fig.js] For 7?^ = 2x10^, the code 
has to be run for 500 ^diff before the dynamo field starts becoming 
'strong' again for the case with ao = 2a c. Due to long compu- 
tational times involved in this exercise we have not continued 
the calculation beyond 60 ^diff • Hence, the determination of sat- 
uration magnetic energy may be inaccurate forRj^ = 2 x 10^. 
Compare this with the case of a simple algebraic quenching of 
the form given in Eq. ^ with g^, = 1. The slopes in the kine- 
matic phase are now almost similar for all /^m within the error in 
the numerical determination of the critical ac. For gc^ = R^, the 
algebraically and dynamically quenched a effects seem to give 
similar dependences on Rj^. It may occur that the two source 
regions may not be spatially separated, so we repeat our cal- 
culations with the a region at r^ = 0.877?© instead of 0.77 Rq 
and obtain the same slope in the relation of the volume averaged 
magnetic energy on Rj^ as in Fig. [5] We also verify from the pro- 
files of field components at two different latitudes, as shown in 
Fig.[6| that the region of strong toroidal field B(p is different from 
the layer where poloidal fields are produced by the a effect. 

For the solutions with dynamical a effect, it may be con- 
cluded from the butterfly diagrams of Fig.[7]that the small-scale 
current helicity au is predominantly negative (positive) in the 
Northern (Southern) hemisphere. The nature of the saturation 
curves of the magnetic energy is strongly governed by the ratio 
of ta and Tcy\. For Rj^ = 20, ta ^ T^cyi = O.SS^diff and so there 
are strong oscillations in the butterfly diagram for au, as shown 
in Fig. ITb, whereas for 7?^ = 200, ta ~ Tcy\ the amplitude of 
oscillations is weak because the au decays at the same rate at 
which it is produced due to the effect of the oscillatory source 
term SB; see Fig.JTp. Similarly for 7?^ = 2 x 10^, the decay 
time ta » T^cyi and so the system of equations is overdamped 
as can be seen from the saturation curve (dashed dotted line) in 
Fig. [2J3 where there are amplitude modulations of the magnetic 
field before it settles to a final saturation value. When the code is 
run longer, we start seeing changes in the parity after t > 40^diff • 
However the magnetic energy and the dynamo period Tcy\ re- 



p. Chatterjee et al.: Magnetic helicity fluxes in spherical sheU dynamos 



10'^ 



10° 



>io-^ 

V 



10-^ 



10" 




Dynamics 
quenching 



{3.)a^(0.72,e)x10\ i?^=20 



10-' 



10" 



10' 



10*= 



10** 



10^ 



10° 



10** 



Fig. 5. Volume averaged magnetic energy scaled with the 
equipartition energy in the saturation phase as a function of 
R^ for dynamical a quenching (triangles + solid) and algebraic 
quenching with g^ = 1 (squares + dashed) and with g^ = Rm 
(cross + dashed-dotted). 




0.65 0,70 0,75 0.80 0.85 0.90 

Fig. 6. Radial profiles of A^ and B^p at two diff'erent latitudes (A) 
in the saturated phase for 7?^ = 2 x 10^. 



main fairly constant even while the system fluctuates between 
symmetric and anti- symmetric parity at an irregular time interval 
(see Fig. [8]). This parity oscillation is absent in the corresponding 
models with algebraic quenching. 

3.2. Secondary dynamo waves 

An interesting result emerges when we repeat our calculations 
with ao = 4a c instead of 2a c for Rj^ = 20. The negative au 
generated in the convectively unstable layer penetrates below 
0.13Rq where aj^ = and drives a secondary dynamo wave 
whose direction of propagation is poleward as compared to the 
primary dynamo wave propagating equatorward. This can be 
seen in the butterfly diagram of B^p at 0.727?© in Fig. [9k. Signature 
of the secondary dynamo can also be seen in the Dutterfly dia- 
gram at 0.8/?©. Even though the secondary dynamo wave is en- 
ergetically powered by the kinematic part of the helical convec- 
tion but the direction of propagation is governed by the sign of 
ckk + <^M. This may be compared with an aQ. dynamo driven 
by a supercritical helicity flux (Vishniac & Cho 2001). This 
mechanism however requires finite initial magnetic field. It may 
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be recalled that we have done calculations with an initial field 
~ 10"^5eq- The diff'erence compared to the case above is that 
the mean field dynamo is not driven by supercritical Vishniac & 
Cho fluxes, but it is governed by a local generation of small-scale 
magnetic helicity. We return to the issue of secondary dynamo 
waves driven by difl'usive magnetic helicity fluxes in Sect. |3.3| 

We also have not observed any evidence of chaotic behaviour 
in the range of magnetic Reynolds number 20 < 7?ni < 2 x 10^ 
for supercritical a < 4a c in agreement with Covas et al. (1998). 
However, if the a efl'ect is highly supercritical, the dynamical 
quenching formula for au is insufficient for dynamo saturation, 
and additional algebraic quenching terms must enter (Kleeorin 
& Rogachevskii 1999). 

3.3. Diffusive magnetic helicity fluxes 

Recently, Brandenburg et al. (2009) showed that catastrophic 
quenching in one-dimensional a^ dynamos can be alleviated by 
introducing a Fickian diff'usive flux in Eq. ([6]) given by 



-KVau- 



(11) 



There was an attempt to calculate the difl'usion coefficient k from 
direct numerical simulations and it was found to be ~ O.Sz/t for 
/?ni ~ 20 (Mitra et al. 2010). For /c = 0, the saturation curves 
in Fig. [2J3 show that the ^sat goes through very low values for 
7?in ~ 2 X 10^ and it takes very long to relax to a steady ampli- 
tude. Next we introduce a diff'usive flux with K(r) = KoT](r) in 



Fig. 10 and obtain ^sat ~ 0.1 ^eq and underdamped behaviour. 



However looking carefully at the corresponding butterfly dia- 
grams (Fig.[TT^,b,c,d) we find a poleward propagating mode due 
to radial diff'usion of the au into the stable layers which other- 
wise was not possible for a very high rjt/rjr ratio. Figures [TT^,f 
show meridional snapshots of sign(B(p)(\B(p\/BQ^)^^^ and au in 
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order to get a clearer idea of the distribution of magnetic fields. 
The poleward propagating mode is now driven by supercritical 
diff'usive helicity fluxes, as opposed to supercritical Vishniac & 
Cho fluxes (see Brandenburg & Subramanian 2005 for examples 
of such behaviour). There exists a /c^ ~ 10"^ for 7?ni = 2 x 10^ 
such that the secondary dynamo fails to operate if atq < Kc and the 
volume averaged magnetic energy decays eventually. It should 
be noted that this threshold for k is highly dependent on R^. For 
instance R^ = 2xl{fi and atq = 10"^ produces a dynamo with 
finite saturation magnetic field and dynamo wave propagation 
governed by au where as for k = 10""^ 77, the dynamo shows 
a runaway growth. An interesting behaviour can be discerned 
from the butterfly diagram of the toroidal field for /?m = 2 x 10^ 
and /Co = 10"^ (Fig 103,c). It appears that the behaviour of the 
dynamo is governed by competition between the poleward prop- 
agating mode and the equatorward propagating mode. The vol- 
ume averaged energy (stars-hline in Fig 10 i) shows correspond- 
ing oscillations long after saturation at an period ~ 5 times the 
period of the equatorward propagating mode. It may be recalled 
that it is well established from direct numerical simulations of 
a^ dynamos that a large-scale magnetic field is easily excited on 
the scale of the system i.e., k~^ for a large kf/ki ratio (Archontis, 
Dorch, Nordlund, 2003). The length scale of the magnetic field 
in Figs 103,c and Figs. 1 1 i,c,e is comparable to k^^, which sug- 
gests that the degree orscale- separation may have become in- 
sufficient to write the electromotive force as a simple multiplica- 
tion, as is done in the expression S = aB - r]J, and that it may 
have become necessary to write it as a convolution, which corre- 
sponds essentially to a low-pass filter (see, e.g., Brandenburg et 
al. 2008). However, we have not pursued this aspect any further. 
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Fig. 9. Butterfly diagrams of the toroidal field (a) and (c) and »„ 
(b) and (d) with a = 4ac for 7?^ = 20. 



3.4. Flux transport Babcock-Leighton Dynamo 

Like in §3.1 we find the critical ckbl required to have a self ex- 
cited dynamo. In this case ac = 5.1 m s"^ for Rj^ = 2x10^. We 
pursue the rest of the calculations with ckbl = 6.0 m s"^ in order 
to avoid producing very large au leading to secondary dynamos 
discussed in §3.1. We should emphasize that Eq. ^ represents a 
first order correction to the a and should be treated with caution 
during its use in supercritical regimes. 

At first we artificially turn off the advective ffux due to 
meridional circulation as well as the diffusive flux only in 
Eq. ([6]), while having them in the induction equations for 5^ 
and A(p. The saturation curve for Rj^ = 2 x 10^ is now over- 
damped whereas the dynamo fails to generate a finite Bsat for 
Rj^ = 2x 10^ even though it initially has the same growth rate. 
On increasing aBi = 10 ms"^ from 6 ms"^ the saturation curve 
for 7?ni = 2 X 10^ also displays overdamped behaviour. This in- 
dicates that the total a in the domain was simply becoming sub- 
critical and the dynamo was not able to sustain itself through the 
saturation phase. We show the distribution of magnetic helicity 
in the meridional plane for in Fig. 13 i, b. Note that aM inside 
the domain is larger for R^ = 2x IcFcompared to 7?ni = 2 x 10^ 
for the same value of ofBL- 

Inclusion of meridional circulation in Eq. ([6]) means that we 
also require a diffusive flux in Eq. ^ to keep the system nu- 
merically stable. A diffusive flux in this equation is known to 
alleviate catastrophic quenching in a^ (Brandenburg et al. 2009) 
as well as aQ dynamos (Guerrero, Chatterjee & Brandenburg 
2010). It is clear from Fig. 12 that the overdamped behaviour 
after the end of the kinematic phase is suppressed due to a dif- 
fusive flux of au which essentially reduces the effective decay 
time for au to much less than Rj^/rjtk^. It may be noted that the 
dependence of the saturation value of the magnetic energy on 
Rjn is now much weaker than the corresponding variation with- 
out fluxes. In presence of diffusive and advective fluxes due to 
meridional circulation in Eq. ^ the small-scale helicity is dis- 
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(a) Interface dynamo with fluxes 
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Fig. 11. Butterfly diagrams of the toroidal field (a) and (c) and 
^m (b) and (d) with a = lac for 7?^ = 2 x 10^ with a /cq = O.Oly/t. 
Meridional snapshots of (e) (8/8^^)^^^ and (f) ajj^x 10^ for the 
same case. 



Fig. 10. (a) Magnetic energy in the domain scaled with the 
equipartition energy for Rj^ = 2 x 10^ indicated in the figure 
for the case of two layered dynamo of Sect. 3.1 with dynamical 
a quenching with a diff'usive flux with kq = 10"^ (star+ solid). 
The same for 7?ni = 2 x 10^ and kq = 10"^ (solid). The satura- 
tion curve for zero fluxes have been shown by the dashed line, 
(b) and (c) show butterfly diagrams for the toroidal field at the 
depths indicated for 7?^ = 2 x 10^ and kq = 10"^. 



tributed through out the convection zone as shown in Fig. 14 i, 
b. It is instructive to compare the weaker magnitudes of au with 



Fig. 13 



The diffxisive fluxes are therefore crucial for operation of a 
successful mean field aQ. dynamo. However an interesting ob- 
servation is the distribution of aM in the concentrated region at 
the lower part of the convection zone (see Fig. 13 ) in contrast to 
Fig. [14] Even though aBi is a surface phenomena, considerable 
magnetic helicity is generated when the meridional circulation 
sinks the poloidal field lines at high latitudes and brings them 
near the tachocline where toroidal fields are generated. 



4. Conclusions 

We have performed calculations for aQ. dynamos in a spherical 
shell for spatially segregated a and Q source regions. The two 
classes of models we have studied resemble the Parker's inter- 
face dynamo and the Babcock-Leighton dynamo. 

In agreement with earlier work, it is not possible to es- 
cape catastrophic quenching by merely separating the regions 
of shear and Of-eff'ect. The saturation value of magnetic energy 
decreases as ~ R^ for both dynamical quenching and the alge- 
braic quenching with ga = Rm for the simple two layer model 
without meridional circulation (Fig. [5]). However we find that a 
richer dynamical behaviour emerges for the cases with dynam- 
ical a eff'ect, in terms of parity fluctuations and appearance of 
'secondary' dynamos (Fig. 8, 9). We do not see evidence for 
chaotic behaviour in the time series of magnetic energy since 
the dynamo period and the saturation energy remains fairly con- 
stant. However this may not be the case in presence of difl'usive 
helicity fluxes which introduce further complexity to the system. 
Addition of difl'usive helicity fluxes relaxes the catastrophic R:^ 
dependence of the saturation magnetic energy (Fig. 10a, 12). An 
interesting 'side-efl'ect' of difl'usive helicity fluxes is the appear- 
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Fig. 12. B^/BI for the flux transport dynamo model of §3.2 for 



Rm = 2x 10^ with /Co = 0.3 (dashed); Rm = 2x 10^ with /cq = 0.3 
(solid); 7?m = 2 X 10^ with /c = (dashed-dotted); 7?^ = 2 x 10^ 
with /Co = (diamond+dashed). 
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Fig. 14. Meridional cross-sections showing the distribution of 
-28.7 toroidal field and aM for a Babcock-Leighton dynamo with MC 
and diff'usive helicity fluxes for R^a = 2 x 10^ at two difl'erent 
epochs. The streamlines of the positive and negative poloidal 
field are shown by solid and dashed lines respectively. 



(b) {B/B )xlO' 
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Fig. 13. Meridional cross-sections showing the distribution of 
toroidal field and aM for a Babcock-Leighton dynamo without 
MC and diff'usive helicity fluxes in Eq. [6]for (a) 7?ni = 2 x 10^ 
and (b) Rj^ = 2x10^ . The streamlines of the positive and negative 
poloidal field are shown by solid and dashed lines respectively. 
Note that the magnetic field has decayed to very small values for 
7?m =2x10^. 



ance of poleward propagating secondary dynamos. However, be- 
cause of the lack of scale separation between the mean field and 
the forcing scale of the helical turbulence we refrain from inter- 
preting this in terms of the poleward migration seen in the Sun. It 
remains to explore the role of the solar wind, coronal mass ejec- 
tions which might help in throwing out the small scale helicity 
from the Sun and thus alleviate catastrophic quenching. The ef- 
fects of Vishniac & Cho fluxes have been investigated and were 



found to be of secondary importance compared to difl'usive he- 
licity fluxes for aQ. mean field dynamos (Guerrero, Chatterjee & 
Brandenburg 2010). 

When both the meridional circulation and the diff'usive helic- 
ity fluxes are artificially shut off' in the helicity evolution equa- 
tion, the dynamo fails to reach significant saturation values, as 
expected (Fig. 12). It is interesting that the Babcock-Leighton 
dynamos, where a is concentrated only in a narrow layer at the 
surface, also produce considerable helicity inside the convec- 
tion zone when the dynamical quenching (Eq. |6]) is employed 
(Fig. 13, 14). 

We have to be cautious about using dynamical quench- 
ing equation for dynamo numbers not very large compared to 
the critical dynamo number. For highly supercritical a, the be- 
haviour of the system begins to be governed by au- We would 
expect that the magnetic field should aff'ect all the turbulent co- 
efficients including both a and rj. However for this analysis we 
have not included an equation for the variation for rj^. This is 
justiffed for the simple two layer model with a lower rjt in the 
region of production of strong toroidal ff elds and a higher 7/t in 
the region of weaker poloidal ffelds. It may also be noted that 
by quenching the diff'usivity inversely with the magnetic energy 
in a nonlinear dynamo model, Tobias (1996) was able to pro- 
duce a bonaffde interface model where the magnetic ffeld was 
restricted to a thin layer at an interface between a layer of shear 
and cyclonic turbulence. However none of the previous interface 
models have used the dynamical quenching equation. 

Unfortunately the direct numerical simulations have not yet 
reached the modest Reynolds numbers used in this paper (~ 10"^) 
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which are still much lower than the astrophysical dynamos. To 
verify if the equation for dynamical quenching works in the same 
way as in a^ dynamos, we need to embark upon systematic com- 
parisons between DNS with shear and convection and mean field 
modelling for aQ dynamos. 
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